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ABSTRACT 

In the present paper we derive the density distribution of dark matter (DM) in a well- 
observed nearby disc galaxy, the Andromeda galaxy. From photometrical and chemical 
evolution models constructed in the first part of the study (Tamm, Tempel & Tenjes 
2008, hereafter Paper I) we can calculate the mass distribution of visible components 
(the bulge, the disc, the stellar halo, the outer diffuse stellar halo). In the dynamical 
model we calculate stellar rotation velocities along the major axis and velocity disper- 
sions along the major, minor and intermediate axes of the galaxy assuming triaxial 
velocity dispersion ellipsoid. Comparing the calculated values with the collected obser- 
vational data, we find the amount of DM, which must be added to reach an agreement 
with the observed rotation and dispersion data. 

We conclude that within the uncertainties, the DM distributions by Moore, Burk- 
ert, Navarro, Frenk & White (NEW) and the Einasto fit with observations nearly 
at all distances. The NEW and Einasto density distributions give the best fit with 
observations. 

The total mass of M31 with the NEW DM distribution is 1.19- IO^^Mq, the ratio 
of the DM mass to the visible mass is fO.O. Eor the Einasto DM distribution, these 
values are 1.28- IO^^Mq and 10.8. The ratio of the DM mass to the visible mass inside 
the Holmberg radius is 1.75 for the NEW and the Einasto distributions. Eor different 
cuspy DM distributions, the virial mass is in a range (6.9-7.9)-10^^Mq and the virial 
radius is ~ 150 kpc. The DM mean densities inside 10 pc for cusped models are 33 and 
16 Mqpc"'^ for the NEW and the Einasto profiles, respectively. Eor the cored Burkert 
profile, this value is 0.06 Mqpc"'^. 

Key words: galaxies: individual: Andromeda, M31 - galaxies: kinematics and dy- 
namics - dark matter. 



1 INTRODUCTION 

Models of hierarchical merging of cold dark matter describe 
rather well the observed general properties of galaxies, large- 
scale structure of clusters, superclusters and their network. 

On the other hand, the observed number of dwarf galax- 
ies seems to be too small in comparison with cold dark mat- 
ter models. But by far the most commonly referred con- 
tradiction is related to disc galaxies. Cosmological simula- 
tions generate dark matter (DM) haloes with central density 
cusps of p ~ or steeper (e.g. Navarro, Frenk & White 
1997; Moore et al. 1999). Observations of dwarf and low sur- 
face brightness disc galaxies have usually shown that shal- 
low central density profiles fit the data better than cuspy 
profiles (Moore 1994; Burkert 1995; Blais-Ouellette, Am- 
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ram & Carignan 2001; Borriello & Salucci 2001; de Blok, 
McGaugh & Rubin 2001; de Blok & Bosma 2002; Salucci, 
Walter & Borriello 2003; Weldrake, de Blok & Walter 2003; 
Simon et al. 2005; Zackrisson et al. 2006; Kassin, de Jong 
& Weiner 2006; Gentile et al. 2004, 2007; Valenzuela et al. 
2007). Observational data are consistent with cuspy density 
profiles only for a few objects (e.g. Swaters et al. 2003; de 
Blok 2005; Simon et al. 2005). 

Discrepancy between the predicted cuspy DM density 
profiles and the observed rotation curves with shallow pro- 
files has been explained with limited resolution of rotation 
curves and non-circular motions of gas and with the state- 
ment that the actual contribution of visible stellar matter is 
poorly known. 

The aim of the present paper is to derive the den- 
sity distribution of dark matter in a nearby disc galaxy, 
the Andromeda galaxy. The galaxy M 31 was selected be- 
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cause (1) photometrical and kinematical (rotation, disper- 
sions) data are measured with sufficiently high resolution in 
order to study the bulge region; (2) velocity dispersions have 
been measured also outside the galactic apparent major axis; 
in addition to stellar kinematics, the kinematics of planetary 
nebulae (PN) is known; (3) direct measurements of metal- 
licities allow to constrain the mass-to-light (M/L) ratios of 
visible matter; (4) independent estimates of the mass distri- 
bution on large scales are available (globular clusters (GC), 
satellites, stream, kinematics of the Milky Way -|- M 31 sys- 
tem). 

In the first part of the study (Paper I) , we constructed 

the photometrical model of M 31 stellar populations on the 
basis of surface brightness profiles in U, B, V, R, I and 
L colours (see also Tenjes, Haud & Einasto 1994). The de- 
rived photometrical model gives us parameters of galactic 
components, colour indices among them. Prom independent 
spectral observations, the metallicity of the stellar content 
is available. The colour indices and the metallicity of each 
component were interpreted with the help of chemical evo- 
lution models to calculate the ages and M/L-ratios of the 
components. In total, the output of the first paper were den- 
sity distribution parameters, ages and M/L-ratios of galactic 
components. 

In the present paper, we apply the results of Paper I 
and construct a mass distribution model of M 31 consistent 
with the measured kinematics. Stellar components and their 

Ad/L's together with the photometrical model give us the 
mass distribution of visible matter. Calculating the rotation 
velocities and velocity dispersions of visible matter with the 
help of the dynamical model (see Tcmpcl & Tenjes 2006) 
the amount of DM can be found, which must be added to 
reach an agreement with the observed rotation and disper- 
sion data. 

In principle, to derive the DM density distribution, it 
is not necessary to take into account stellar velocity dis- 
persions. It is sufficient to compare the gas rotation curve 
with the calculated derivatives of the gravitational potential. 
This simple and straightforward method is frequently used. 
Of course, before calculations the gas rotation curve must be 
corrected for possible expansion velocities, non-circular mo- 
tions and dispersions (see Gentile et al. 2007). However, one 
of the most uncertain aspects of this method is the M/L of 
the visible matter. Chemical evolution models involve sev- 
eral insufficiently constrained parameters. It is possible to 
decrease degeneracy by using 3 5 different colour indices, 
but due to measurement uncertainties, colour indices are of- 
ten controversial. The comparison of the results with addi- 
tional and independent observations (stellar rotation curve, 
velocity dispersions along several slit-positions) allows us to 
constrain the distributions of visible matter and thereafter 
DM. 

In the present model, we assume that a galaxy is a su- 
perposition of ellipsoids of rotation with different flattening. 
We assume that stellar populations have both rotation and 
dispersion components. The velocity dispersion ellipsoid of 
visible matter is triaxial. DM has spherical density distribu- 
tion and is collision-free. 

The general scheme of Paper I and Paper II (present 
paper) is given in Fig. 1. The photometrical model and the 
chemical evolution model constitute the first part of the 
study and depend very strongly on each other. The output of 
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Figure 1. The principle scheme of the construction of the M31 
model in the framework of Paper I and Paper II (the present 
paper). The arrows show dependencies between different model 
construction stages. The dashed arrows represent weaJcer impact. 



the first part (component parameters and M/L-ratios) is the 
input to the second part of the study. At the second stage, 
we use the kinematical model and narrow down the output 
of the first part. Finally, we get a self-consistent model and 
derive the DM density distribution, using a possibly wide 
variety of observational data. 

In Section 2, we describe the observational data which 
we have used. Kinematical data consist of gas rotation, stel- 
lar rotation and dispersions, rotation velocities and disper- 
sions of PN and velocity dispersions of individual red gi- 
ant branch (RGB) stars. In Section 3, the dynamical model 
and an algorithm for calculating the DM distribution are 
described. In Sections 4 and 5, visible and DM density dis- 
tributions are derived and compared with available observar 
tional data. In Section 6, the results are discussed. 

The general parameters of M 31 used in the present pa- 
per are: the system velocity is -300kms~^ (de Vaucouleurs 
et al. 1991; van den Bergh 2000), the galaxy inclination angle 
is 77.5° (Walterbos & Kennicutt 1988; de Vaucouleurs et al. 
1991), the major axis position angle is 38.1° (Walterbos & 
Kennicutt 1987; Ferguson ct al. 2002) and the distance has 
been taken 785 kpc (McConnachie et al. 2005), correspond- 
ing to the scale 1 arcmin = 228 pc. Throughout the paper 
we use the ACDM cosmological model with = 0.3 and 
Ho = 75kms"^Mpc-\ 



2 OBSERVATIONAL DATA 

In Paper I we described photometrical and chemical com- 
position data used to construct the photometrical model. 
In the present section, we describe kinematical data, used 
for the construction of the dynamical model. From a large 
variety of kinematical information, we present only the ob- 
servations, that can be compared with our model output. 
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Figure 2. Circular velocity, derived from Hi and Hll observa- 
tions. The filled circles are averaged rotational velocities obtained 
using the data from earlier observations (faint points). The filled 
squares are recent Hi observations by Carignan et al. (2006). See 
the text for more information. 



Table 1. References for rotation and dispersion data along the 
major and the minor axis. 



Reference 


-Rmax [kpc] 


Zmax [kpc] 


McElroy (1983) 


2.27 


1.98 


Kormendy (1988) 


0.185 




van der Marel et al. (1994) 


0.234 


0.16 


Kormendy & Bender (1999) 


0.18 





2.1 Hi and Hll observations 

We constructed the gas rotation curve of M31 on the basis 
of H II observations by Rubin & Ford (1970) and H I observa- 
tions by Kent (1989) and Braun (1991). Observational data 
points from these studies were averaged, and the resulting 
rotation velocities are presented by filled circles in Fig. 2. To 
calculate the error-bars, a bootstrap method is used, com- 
bining the bootstrap errors with observational errors. 

In outer parts of M31, an extended Hi rotation curve 
has recently been measured by Carignan et al. (2006) . They 
measured the rotation curve outside 90 arcmin, between 20 
and 35 kpc, and recomputed also rotation velocities, using 
earlier Hi data by Unwin (1983). Rotational velocities by 
Carignan et al. (2006) are showen by filled squares in Fig. 2. 
The derived rotation curve represents gas kinematics from 
1-35 kpc. 



2.2 Stellar velocity and dispersion observations 

The intrinsic structure of the nucleus of M31 is not clear - 
it may have a ringlike and/or a double or even a triple struc- 
ture (Lauer et al. 1993; Gerssen, Kuijken & Merrifield 1995; 
Peiris & Tremaine 2003; Salow & Statler 2004; Bender et al. 
2005). In The context of the present paper, the contribution 
of the nucleus is unimportant, we look at regions outside 
0.01 kpc and we exclude the nucleus from our model. 

Stellar rotation velocities and velocity dispersions in the 
central regions of M 31 have been measured by several au- 
thors (see references in Tenjes et al. 1994). To model the con- 
tribution of visible and DM to the total mass distribution, 
it is useful to possess kinematical data outside the galactic 
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Figure 3. Stellar dispersion observations. The faint points cor- 
respond to observations from different authors. The filled circles 
denote the averaged profile. 
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Figure 4. Stellar rotational velocities. The faint points are ob- 
servations from different authors. The filled circles denote the 
averaged profile. 



apparent major axis and at largest distance intervals (see 
Tempel & Tenjes 2006). References to stellar kinematical 
data used in the present paper are presented in Table 1 
which also gives maximal distances along the major and 
the minor axis. McElroy (1983) measured stellar rotation 
and dispersions up to the distances of 2.3 kpc from the cen- 
tre along 12 cuts (at every 15°) that penetrate the galaxy 
centre. None of them was exactly along the galactic major 
(minor) axis, and for this reason, we averaged two cuts at 
6.8° and 8.1° from the major (minor) axis. Observations by 
other authors were sufficiently close to the major (minor) 
axis. 

Velocity dispersions and rotation velocities along the 
major axis are presented in Figs. 3 and 4, respectively. 
Individual observations have been averaged; for error-bars 
a bootstrap method was used, combining bootstrap errors 
with observational errors. Outside the major and the minor 
axis, the observations by McElroy (1983) were combined into 
four cuts at 23°, 38°, 53° and 68° from the major axis and 
they are shown by filled circles in Fig. 10. 



2.3 Planetary nebulae observations 

Kinematical observations for PN have been published by 
Halliday et al. (2006); Merrett et al. (2006). Halliday et al. 
(2006) used 50 nebulae for every data point and Merrett 
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ct al. (2006) used roughly 40 nebulae. For every data point, 
they have found the averaged velocity and dispersion near 
the galactic major axis. 

In our model, we handle the kinematics of PN in a sim- 
ilar way to the kinematics of the stellar component. As a re- 
sult, the stellar observations give us the kinematics in inner 
parts of the galaxy and PN observations give the kinematics 
mainly in outer parts. In overlapping regions both data are 
in accordance. 



2.4 Observations of red giant branch (RGB) stars. 

Reitzel & Guhathakurta (2002) observed the kinematics of 
RGB stars along the minor ax:is. At a distance of 19 kpc from 
the centre, the calculated velocity dispersion was found to 
be 150t^E]kms"^ Kalirai ot al. (2006) measured 106 RGB 
stars at 12 kpc along the minor axis and the derived disper- 
sion was 90 kms~^. The most recent observations by Gilbert 
et al. (2007) involved distances in the range of 9-30 kpc along 
the minor axis. They measured RGB stars in different fields 
and the averaged dispersion in this region is 129kms~^. The 
results from different fields and errors are presented in Fig. 9. 

Chapman ct al. (2006) made observations of RGB stars 
along the major axis in the halo region (10-70 kpc). They 
made spectroscopic observations for 827 stars. The main 
conclusion was that the halo sample shows no significant 
evidence for rotation. They mcEisured the central dispersion 
152kms~^, which decreases with radius. The observed dis- 
persions by Chapman et al. (2006) are higher than stellar 
or PN observations. There is no strict explanation for that, 
but probably the linc-of-siglit velocity component is not cor- 
rectly subtracted, and therefore the dispersions are higher. 
Also, there is a possibility that RGB stars have systemati- 
cally higher dispersions. Unfortunately, in our final model we 
cannot compare these observations with our model output. 



3 A DYNAMICAL MODEL 

3.1 Calculation of stellar velocities and 
dispersions 

Details of the dynamical model described below can be 
found in Tempel & Tenjes (2006). However, to improve the 
readability of the paper, a brief outline is given here. 



3.1.1 Basic formulae 

Let (R,z,6) be cylindrical coordinates and a = 
^JB? + z?/(p , where q is the axial ratio of isodensity el- 
lipsoids. Knowing spatial luminosity densities of the compo- 
nents li{a) and ascribing a M/L-ratio to each component /, 
{i indexes the bulge, the disc, the stellar halo and the outer 
diffuse stellar halo; see Paper I), we obtain the spatial mass 
density distribution of a galaxy 



p{a) = ^ fili{a) + pDM(a) 



(1) 



(pDM(ffl) is the DM density). On the basis of spatial mass 
density distributions, derivatives of the gravitational poten- 
tial 1^ and |j can be calculated. 



In stationary coUisionlcss stellar systems with axial- 
symmetry the Jeans equations in cylindrical coordinates can 
be written in a convenient form for further calculations 



dR 

dz 



fl — ke 9k\ 



ozj az 



+ 



dR 



9$ 
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where 

1 - 



tan(2a) 



k =^ 



ke^4- 



(3) 



(4) 



(5) 



a is the angle between the major axis of the velocity disper- 
sion ellipsoid and the galactic plane. For each component, 
the rotation velocity has been taken Vg = pVc, where K is 
circular velocity and /3 is a constant specific for each sub- 
system. 

The Jeans equations (2) and (3) contain unknown func- 
tions kz, kg and a. In our model, the phase-density of a 
system is a function of three integrals of motion. Under this 
assumption, wc can derive the unknown functions. 

From the Jeans equations (2) and (3) the dispersions 
along the coordinate axis (an, az and ag) can be calculated 



oo r 

paUR, «)=(!- /?') I P^^^ exp I p{r\ z)dr 



2 , d^{R,z') 

paz{R,z) = I p- 



Qz 



exp ; g{R, z*)dz' 



dz' , 



2 1 2 2 kg 2 

ag = kgaji or ag = -r-a^, 



where 

i-kg , a 
p = — - — h 



9 



dr, (6) 

(7) 
(8) 

R ' dz'- ^ ~ R ' dR' 
The details of derivation and some restrictions can be 
found in Tempel & Tenjes (2006). 

3.1.2 Line-of-sight dispersions 

Calculated velocity dispersions cr^, a^ and ag cannot be 
directly compared with observations. The calculated disper- 
sions must be projected to the linc-of-sight. Designating Q 
as the angle between the line-of-sight and the galactic plane 
(see Fig. 5), the line-of-sight dispersions af^^ are 

2 2 2 , 2-2 /-\ 

cTios = cr» cos + a^ sin B, 



(10) 



where 



cr* = ag— + an I 1 - ^ 



R^ 



i?2 



(11) 



The next step is to integrate individual line-of-sight ve- 
locity components along the whole line-of-sight 



al,{X,Y)= ^ 



L{X,Y) 



Y: [KR,z,)aUR,z.)] 



R- 



coseV^2 _X2 



-dR, (12) 
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Figure 5. Panel (a) gives the angle between the line-of-sight and 
the galactic plane. Panel (b) shows our galactic coordinates X 
and Y in the sky. 



Zl,2 = 



Y 
sin O 



± ^/R? -xA tan e 



(13) 



where 1{R, z) denotes galactic spatial luminosity density, and 
L{X, Y) is the surface luminosity density profile. 

Equation (12) gives the line-of-sight dispersion for one 
galactic component. Our model consists of several compo- 
nents and wo must sum over all components considering the 
surface luminosity profile 



Y,{U{X,Y) [aL(X,F)].} 



al^,{X,Y) = 



(14) 



where i denotes the subsystem and the summation is taken 
over all subsystems. 

3.1.3 Line-of-sight stellar velocities 

In an axisymmetric system, wc have only one velocity com- 
ponent, the rotational velocity Ve. According to the defini- 
tion of the circular velocity, we can write for each component 
(Binney & Tremaine 1987, equation (2.19)) 



dR 



(15) 



where = + z'^x^{I^ — e'^x'^)~^ and e = ^J\ — is 
eccentricity. The rotational velocity cannot be compared di- 
rectly to the observed velocity. We use the same projection 
method for stellar velocities as we used for dispersions. The 
projected line-of-sight velocity is 



Vios = Ve—cosQ. 
H 

The observable velocity Vobs is 

Y.{Li{x,Y) [M„t(x,y)].} 



(17) 



where Vint is the same as the right-hand side of equa- 
tion (12), only cTiQg = Vios. The summation is taken over 
all subsystems. 

3.2 Calculation of dark matter distribution 

Total gravitational potential is a sum of the potentials 
caused by visible matter and DM. In terms of circular ve- 
locities, 

V^°' = V^''' + V^y,. (18) 

Photometrical and chemical evolution models give us 
the spatial mass density of visible matter (see equation (1)) 



and thereafter I4 If g£is velocity dispersions are small 
in comparison with rotation velocities and we neglect also 
non-circular motions, gas rotation velocities can be identified 
with total circular velocities and equation (18) gives us the 
spatial mass density distribution of DM 

and inner mass of DM 

2 



Mjym{R) = 



G 



(20) 



where G is the gravitational constant. 

Unfortunately, this straightforward calculation is too 
simplified. In the most interesting bulge region, ignoring of 
non-circular motions and dispersions of the gas is not jus- 
tified. For this reason, total circular velocities remain un- 
known and wc must use another way to calculate the DM 
density. 

In a simple form, the first Jeans equation (2) can be 
written as 



V! = Vi + alf{R), 



(21) 



where f{R) is a specific function. Applying (21) for total 
matter (visible -|- DM) , we have 



V^°' = Ve^°' +aTf{R). (22) 
Substituting the last equation into equation (18), we get 



hot 



(23) 



where (Jfl.°* and Vg°'^ arc the observed velocity dispersion and 
the observed rotational velocity, respectively. The circular 
velocity of the visible matter results from the component 
parameters derived in Paper I. It is reasonable to assume 
that the function f{R) is the same for total matter and 
for visible matter (DM is collision- free) . Taking this into 
account, 



A-is 



K"''= = + (24) 

Expressing from the previous equation f{R), we can write 



(16) f{R) 



(25) 



2 . 

where V^^^ and 0-?'= are the visible matter rotational veloc- 
ity and dispersion, respectively. We can find those quanti- 
ties, using the method described in Section 3.1. Substituting 
equation (25) into equation (23), we finally get the DM cir- 
cular velocity 

Vdm = Vp°' - V^"" + 4? [V^'" - V^'^) ■ (26) 

Using equations (19), (20) and (26), we find the DM density 
distribution and the inner mass of DM, taking into account 
the rotation and dispersion data. 

Rotational velocities and dispersions in equation (26) 
must be determined in the galactic equatorial plane as a 
function of radius. 

The Jeans equation in form (21) is not entirely correct, 
because we did not consider dispersion derivatives. To have 
the correct form, an additional term must be added to the 
right-hand side 
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-R 



9p 
9i? 



an 



(27) 



Taking this correction into account, also a correction to the 
right-hand side of equation (26) must be added 



R 



dR dR 



+ RcTob 



1 9pvis 1 9/5t. 



dR ptot 3-R 



(28) 



This additional term is important only on a small scale. 
We cannot adequately consider the derivatives of dispersions 
as numerical differentiation amplifies noises in observational 
data too much, and therefore we do not use this term in our 
calculation. 



4 MASS DISTRIBUTION MODEL OF M 31 

We assume M 31 to consist of the following subsystems: a 
bulge, a disc, a stellar halo, an outer diffuse stellar halo and 
a DM halo. Reasoning of the component selection (except 
DM halo) and their parameters resulting from the photo- 
metrical model are given in table 1 in Paper I. Luminosity 
distributions of the components are transformed to mass 
distributions on the basis of AI/L's resulting from chemical 
evolution models. Chemical evolution models give the com- 
ponents most probable AI/L in different colours and possible 
M/L ranges. 

As most kinematical observations have been made at 
wavelengths close to the V-band, the ^/-luminosities are 
used in weighting the components if superposition is needed 
(Equations (14) and (17)). Component parameters most rel- 
evant in context of the present paper are given in Table 2. 
According to Nieten et al. (2006), we add to the disc a total 
neutral gas mass 5.2 • IO^Mq. 

Dynamical models can further constrain the permitted 
mass intervals. In our dynamical modelling process, we tried 
many different models, and, in most cases, the main con- 
flict was that chemical evolution models predict larger bulge 
masses than dynamical models permit. The masses of other 
components remain between the permitted values. In our fi- 
nal dynamical model, masses for the disc, the stellar halo 
and the outer diffuse stellar halo are in accordance with the 
masses from the chemical evolution model. To achieve the 
correct concentration parameter for the NFW dark matter 
halo, we slightly decreased the best disc mass resulting from 
the chemical evolution models. The bulge mass in our fi- 
nal model is also slightly smaller than the minimal stellar 
mass predicted by the chemical evolution model. The smaller 
bulge mass was needed for an agreement with the observed 
velocity dispersion data in the bulge region. 

In addition to masses, dynamical models contain a pa- 
rameter /3. This parameter was, in fact, determined for each 
component on the basis of fitting the stellar and PN rota- 
tion velocities with the model. On the other hand, an in- 
creasing of rotation velocities causes a decrease of velocity 
dispersions. In this way also dispersions influence the de- 
termination of f3. Final values of f3 parameters are given in 
Table 2. 

Circular velocities for visible and DM components in 
the final model are given in Fig. 6 as dashed lines. For com- 
parison, observed gas rotation velocities are given by filled 
circles and filled squares. Our model is in good agreement 
with observations. 
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Figure 6. Gas rotational velocity. The filled circles are an aver- 
aged profile (see Fig. 2), the filled squares show the rotation curve 
from Carignan et al. (2006). The solid line denotes circular veloc- 
ity derived from our model. The dashed lines show the rotational 
velocity for different components. 
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Figure 7. PN and stellar velocity along the major axis. The 
filled triangles and the diamonds are PN rotational velocities, the 
filled circles are averaged stellar velocities. The solid line is our 
calculated stellar velocity. 
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Figure 8. Dispersions along the major axis. The filled triangles 
and the diamonds are PN dispersion observations. The filled cir- 
cles are averaged stellar dispersions. The solid line is the disper- 
sion from our model. 
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Table 2. Calculated model parameters for stellar components. 



Popul. 



ao 



N 



M/Lv 



M 



Bulge 


1.45 


0.64 


0.6 


4.2 


.2+1-0 


3.39 


0.25 


Disc 


2.02 


9.3 


0.05 


0.7 


o ,+1.2 
•^■1^-0.7 


5.57 


0.98 


Halo 


0.61 


4.0 


0.5 


2.7 


2 

^■^-0.7 


1.76 
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Masses and luminosities are in units of 101°Mq and IO^^I^q re- 
spectively; component radii are in kpc; M/Ly is expected M/L- 
ratio from chemical evolution models and M is mass used in our 
dynamical model. 
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Figure 9. Stellar dispersions along the minor axis. The filled cir- 
cles and the diamonds are stellar dispersions. The pentagon, the 
square and the triangles are individual RGB stars observations. 
The solid line is our modelled dispersion. 



Stellar rotational velocity is given in Fig. 7. In the in- 
ner region [R < 0.02 kpc), the observational velocities are 
higher than the calculated velocities. M31 has a fast rotat- 
ing nucleus and the increase of velocities can be explained 
by using the nucleus. Our model does not include the nu- 
cleus and therefore we have a lower velocity profile in this 
region. The averaged observed stellar rotational velocities 
in the outer region {R > 0.7 kpc) are smaller than the cal- 
culated velocities. We do not have a good explanation for 
that. However, the observations are very uncertain in this 
region (see Fig. 4) and some observations are in agreement 
with our calculated results. We have a satisfactory fit with 
PN observations. A small disagreement is only in the outer 
regions, where the calculated velocities are slightly smaller 
than the observed ones. 

Calculated velocity dispersions along the major axis are 
given in Fig. 8. It is seen that we have a good fit with ob- 
servations. Only in the outer region {R > 5 kpc, where only 
PN observations are available) the calculated dispersions are 
slightly smaller than the observed dispersions. The proba- 
ble reason for that is that PN observations do not lie exactly 
along the galactic major axis. The second reason is that the 
observed dispersions include a component caused by rota- 
tion at different positions along the line of sight, which in- 
creases the observed dispersions by up to twenty kilometers 
per second. 

Calculated velocity dispersions along the minor axis 
are presented in Fig. 9. In addition to usual stellar dis- 
persion measurements, we have plotted in Fig. 9 also in- 
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Figure 10. Dispersions along slits between the major and the 
minor axis. The filled circles are dispersions from McElroy (1983). 
The solid lines are our modelled dispersions. The distance r is 
along the slit. 



dependent dispersion measurements for RGB stars by Re- 
itzel & Guhathakurta (2002); Kalirai et al. (2006); Gilbert 
et al. (2007), extending to larger distances z from the galac- 
tic plane. The dispersions calculated from our best-fitting 
model are in good agreement with observations of RGB 
stars. In the intermediate regions, the calculated dispersions 
lie higher than the observed ones. It is possible to have a bet- 
ter fit by using the elliptical DM halo, but still the calculated 
dispersions are greater than the observed ones. 

In Fig. 10 the observed and calculated dispersions along 
the lines between the major and the minor axis (at 23° , 38° , 
53° and 68° from major axis) are presented. The observa- 
tions are with high scatter and they are not symmetric with 
respect to the galactic centre and are thus uncertain. 



5 DARK MATTER DISTRIBUTION IN M 31 

One of the aims of Paper I and the present paper is to de- 
rive constraints to the possible DM density distribution and 
to compare these constraints with several known analytical 
spherical DM density profiles. 

First, we chose for comparison a typical cored density 
distribution (Burkert 1995) 

Po 



PBurkcrt('") 



(29) 



where po is central density and Vc is characteristic radius 
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The cored isothermal density distribution profile is rather 
similar to the Burkert profile (with the exception of the outer 
parts) and we left it aside. Next we also used different cuspy 
density distributions: the Moore profile (Moore et al. 1999) 



Table 3. Dark matter profiles. 



PMoore(f) 



and the NFW profile (Navarro et al. 1997) 



PNFw(f) 



Pc 



(30) 



(31) 



In the Moore and NFW profiles, pc is a density scale param- 
eter. 

Later Navarro et al. (2004) have found instead of the 
NFW profile another density distribution law which fits over 

a wider distance range. According to Morritt et al. (2006), 
we refer to this equation as "Einasto's r^^" model" (Einasto 
1965, 1968, 1969). Following Merritt et al. (2006), we replace 
the exponent a in Navarro et al. (2004) equation by 1/n 



PEinasto = Pc exp < —dn 



l/n 



1 



(32) 



where n is in principle a free parameter. According to N- 
body simulations, we can take n = 6.0 ± 1.1 (Navarro et al. 
2004). The term d„ is a function of n in a way that Pc is 
the density at rc defining a volume containing half of the 
total mass. The value of dn can be well approximated by 
the expression (Merritt et al. 2006) 



d„ « 3n- 1/3 + 0.0079/n. 



(33) 



The NFW density profile has extensively been studied 
and a useful characteristic for this profile is the concentra- 
tion parameter, cvir, defined as the ratio between the virial 
and inner radii. 



Rvir/r 



(34) 



The outer, virial radius Jivir of a halo of the virial mass Mvir, 
is defined as the radius within which the mean density is Avir 
times the critical density of the universe (pcrit = SHq /8itG): 



47r . 

Mvir 



r Pcrit -^^ii 



(35) 



In the ACDM cosmological model, the local value is Avir — 
337. 

Using N-body simulations, the correlation between the 
halo virial mass and the concentration parameter for the 
ACDM model has been found by Bullock et al. (2001); Wech- 
sler et al. (2002). For an Andromeda-sized halo, the concen- 
tration parameter is in the range 10 < Cvir < 20. In our 
model, the concentration parameter is 11.7, which is in ac- 
cordance with N-body simulations. In principle, in the case 
of Moore's law, the concentration index should also be con- 
sidered. Unfortunately, the dependence of the concentration 
index of Moore's law on DM mass is unknown and we cannot 
take it into account. For the Einasto law, the consentration 
index in our model is siniiliar to that for NFW law. 

Following the scheme presented in Section 3.2, we can 
estimate the allowed region for the DM density distribution 
in M31. 

In Fig. 11 the derived DM density distribution as a func- 
tion of radius is given. To achieve the maximum range in 
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Parameter n has been taJsen 5.8. 



radius, we used three different types of data: stellar obser- 
vations, PN observations and the gas rotation curve. Prom 
kinematical data (Fig. 7 and 8) we see that PN give us the 
mass distribution in an intermediate range 0.2 20 kpc. Using 
stellar observations, we can obtain the mass distribution in 
inner region. In the outer region (up to 35 kpc), we estimate 
the mass distribution using the gas rotation curve. Using 
three different types of observations we obtain the DM mass 
distribution in the range 0.02 35 kpc. The mass distributions 
derived from PN data have an intersection with stellar and 
gas data. In interlapping regions all three estimates are in 
accordance. The error-bars (filled areas) in Figs. 11 and 12 
are derived using kinematical observational errors only. The 
real error-bars arc slightly larger, because we must add er- 
rors from photometrical observations and intrinsic errors of 
the model. 

Derived in this way density distribution has been fitted 
with the analytical profiles referred above. Best-fitting pro- 
files are plotted also in Fig. 11 and 12, the corresponding 
parameters of the profiles are given in Table 3. In addition, 
virial masses Mvir and virial radii iivir for each profile are 
given. In our final model (described in the previous section), 
we have used the Einasto profile. 

On the basis of Fig. 11, the inner mass of the DM com- 
ponent was calculated and plotted in Fig. 12 as a function of 
radius. A number of researchers have used different objects 
to estimate the total mass inside a given radius or the total 
mass of M 31. In Table 4 we present some mass estimates at 
different radii with error-bars. According to Geehan et al. 
(2006), most realiable estimates are from Evans & Wilkin- 
son (2000), using the PN and the GC data and an estimate 
by Geehan et al. (2006) on the basis of satellite galaxies. At 
a distance rmax = 125 kpc three independent mass estimates 
have been made with the mean value 7.5-lO"M0. The re- 
ferred mass estimates in Fig. 12 correspond to dark matter 
masses: from total masses given in Table 4 the stellar mass 
was subtracted. 

When comparing the derived DM distribution with 
most common analytical DM density distributions (Figs. 11 
and 12), it is seen that within uncertainties nearly all distri- 
butions fit with the limits from observations. On the other 
hand, it is also seen that the Moore profile lies in central re- 
gions slightly outside the upper limits of the uncertainties; 
the Burkert profiles lie rather close to the lower limits of the 
uncertainties. The Burkert fornmla predicts also a slightly 
smaller total and virial mass than the cuspy profiles. It is 
seen that the NFW and the Einasto density distributions 
give the best fit. 
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Table 4. Mass estimates for M31 (out to a radius Tmax)- 
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Figure 11. DM density profile. The filled areas are density pro- 
files derived form stars, PN and gas data. Different lines represent 
different DM profiles. In our final model we have used the Einasto 
profile. 
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Figure 12. Integrated DM density profile. The filled areas and 
the lines are the same as in Fig. 11. The points are mass estimates 
from different authors using different objects. 



6 DISCUSSION 

In Figs. 11 and 12 we tried to estimate all possible uncer- 
tainties. The contribution of errors in photometrical surface 
brightness profiles is small. There is some uncertainty in 
component subtraction but its infiuence due to integration 
along the line-of-sight is relatively small. The most signifi- 
cant uncertainty arises from the stage of using chemical evo- 
lution models. First, different colour indices are not in good 
mutual agreement. Second, chemical evolution models in- 
volve several free parameters (initial mass function (IMF), 
age, star formation duration, contribution of stellar winds 
etc.). However, it was possible to decrease the influence of 
all these uncertainties in our final results by using several 
colour indices {U - B), {B -V), {V - R), (R-I) and (I-L). 
In addition, we used independent metallicity measurements 
allowing us to significantly narrow the possible models. The 
uncertainties resulting from chemical evolution models were 
estimated by constructing a large amount of different models 
with different initial conditions. All this allows us to ascribe 
errors to the final AI/L and thereafter to the M values of 
visible matter for all components (see Table 2). 

We derived that within uncertainties nearly all the com- 
monly used analytical DM density profiles fit with limits 
from observations at nearly all distances. The best-fitting 
models, however, are the NFW and the Einasto density pro- 
files. Usually in similar studies it is concluded that the NFW 
profile does not fit with observations. 

A common approach in this kind of studies at present 
is to select galaxies without the bulge component or with 
a small bulge (see e.g. Simon et al. 2005; Kassin, de Jong 
& Weiner 2006; Gentile et al. 2004, 2007). The gas rotation 
curve is corrected for non-circular motions and expansion 
as well as possible. Colour indices of the disc are derived 
from surface brightness distributions. Using a colour-M/L 
relation (usually by Bell & de Jong (2001), updates Bell et al. 
(2003)) the mass distribution of visible matter is derived. 
Identifying gas rotation velocities with circular velocities it 
is possible to calculate the DM density distribution. 

To derive the M/L of visible matter, quite frequently 
only one or two colour indices are used. This may involve 
significant uncertainty in the density distribution of visible 
matter and thereafter in DM. 

When modelling SO-Sa-Sbc disc galaxies with a signif- 
icant bulge component, it is necessary to take into account 
a detailed luminosity distribution of the bulge and the disc, 
finite ellipticities of these components, and both the rota- 
tion and the dispersion components in stellar kinematics. 
When ignoring stellar velocity dispersions, the calculated 
gravitational potential is smaller. Dispersion components in 
kinematics of a typical spiral galaxy are more dominating 
in the central parts of the galaxies (see Fig. 13), and thus 
in the central parts the gravitational potential is underesti- 
mated and the resulting DM density profile is fiatter near 
the centre. 

Even in a sample of low mass disc, the power-law index 
in the central regions of galaxies varies highly (Simon et al. 
2005). Thus it is worth not to limit DM density studies with 
bulgeless low surface brightness (LSB) galaxies but to model 
also galaxies of intermediate and high luminosities and of 
different morphological types. In LSB and late morpholog- 
ical type galaxies star formation is delayed. This may be 
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Figure 13. Gas rotation, stellar rotation and stellar dispersion 
components of a typical spiral galaxy as a function of the galactic 
radius. 

caused by insufficiently deep potential well of incompletely 
formed DM haloes. Thus, it is not surprising that DM den- 
sity profiles in these galaxies are dominantly not cuspy. In 
addition, recent study by Pizzella et al. (2008) rises the ques- 
tion about the reliability of the use LSB galaxies as a tracer 
of mass distribution, in particular in the central regions of 
LSB galaxies. 

Modelled dispersions along the minor axis (Fig. 9) and 
between the major and the minor axis (Fig. 10) do not have a 
very good fit with the observed stellar dispersions. Unfortu- 
nately, this bulge region has been observed only by (McElroy 
1983) and there is a possibility that these observations in- 
clude some systematic errors that have not been taken into 
account. On the other hand, our models have several free 
parameters. In our final model, we use spherical DM den- 
sity distributions. We have tried elliptical DM profiles and 
in some aspects the elliptical DM profile gives slightly bet- 
ter results. In our model we have also parameter B, which 
determines the rotational velocity for stellar components. In 
our model, /3 is a constant for each subsystem, but princi- 
pally might be a function of R and z. As we do not have 
good quality stellar rotation curves outside the galactic ma- 
jor axis, wo do not know how l3 varies. The last but not 
least uncertainty in our model is the orientation of the tri- 
axial velocity ellipsoid. In our model, the velocity ellipsoid 
is radially elongated and lies under the angles ^ 30° with 
respect to the galactic equatorial plane (Tempel & Tenjes 
2006). The velocity ellipsoid orientation is quite uncertain 
and the orientation can affect the modelled dispersions quite 
significantly. In one of our following papers we analyse this 
kind of dynamical models in greater detail. 

Our derived total mass inside 100 kpc (5.6- lO^M© for 
NFW and 6.2 • 10"Mq for the Einasto model) is compa- 
rable with mass estimates from Klypin, Zhao & Somerville 
(2002) (8.5 • 1O"M0), Geehan et al. (2006) (6.0 • lO^M©) 
and Seigar, Earth & Bullock (2006) (6.5 ■ 10"Mq). The 
ratio of visible matter to total matter (inside 100 kpc) in 
our model (0.10 for the NFW and 0.092 for the Einasto 
model) is in good agreement with the results from Klypin 
et al. (2002) (0.105). On the contrary, the ratios from Gee- 
han ct al. (2006) (0.173) and Seigar ct al. (2006) (0.152) are 
somewhat higher than ours: the main difference is that they 
have almost twice as massive a disc as in our model. Gee- 
han et al. (2006) admit that they have degeneracy between 



the disc and dark halo components: reducing the disc mass 
and increasing the dark halo mass we get better fit with 
our model. The ratio of the DM mass to the visible mass 
inside the Holmberg radius (~ 25 kpc for M31) is 1.75 for 
the NFW and the Einasto distributions. 

The total mass of M 31 with the NFW DM distribution 
is 1.19 • IO^^Mq, the ratio of the DM mass to the visible 
mass is 10.0. For the Einasto DM distribution these values 
are 1.28-10^'^ Mq and 10.8. Total masses are calculated inside 
400 kpc, which is half way from M 31 to the Milky Way. The 
total mass calculated from our model is in good agreement 
with the total mass calculated by Evans & Wilkinson (2000) 
on the basis of M 31 satellites. 

On the basis of the proper motions of the Local Group 
galaxies IC 10 and M 33 and assuming that these galaxies 
are bound to M31, Brunthaler et al. (2007) calculated a 
lower limit for the mass of M31 of 7.5 • IO^^Mq. This mass 
estimate is in accordance with our mass estimate based on 
cuspy DM profiles. The cored Burkert profile in our model 
gives slightly smaller mass than this lower limit. 

In central regions the density distribution of the DM can 
be approximated as p ~ r~^. From Fig. 11 we can conclude 
that 0^/3^ 1.5. On the basis of a solution of the .leans 
equations Hansen (2004) derived that 1^/3^3. Resulting 
common region is rahter narrow 1 ^ /3 1.5. 

For local dwarf spheroidal galaxies Gilmore et al. (2007) 
derived the central characteristic DM density (averaged 
over a volume of radius 10 pc) within a cuspod model 
-6OM0pc"^. For our models, the corresponding values are 
33 for the NFW profile and 16 for the Einasto profile. For the 
cored Burkert profile, this value is 0.057 Mqpc~^. It is inter- 
esting to note that the central characteristic density of the 
DM for cuspy profiles is nearly independent of the galactic 
mass. This value is close to the upper limit of the central DM 
density range of nine local galaxies (l-4)-10~^* gcm~^ = 
(0.015-0.059) M0pc"^ derived within the cored DM mod- 
els by Borriello & Salucci (2001). For four galaxies at mean 
redshifts (z) ~ 0.9, Tamm & Tenjes (2005) derived the cen- 
tral DM density (0.012-0.028) Mqpc"^. This lower value at 
higher redshift may hint to the cosmological evolution of 
DM densities. 



7 CONCLUSIONS 

In the present paper we have derived the density distribu- 
tion of DM in a most well-observed nearby disc galaxy, the 
Andromeda galaxy. In the first part of the study (Paper I) , 
we construct the photometrical model of M 31 stellar pop- 
ulations on the basis of surface brightness profiles in U , B, 
V , R, I and L colours. The derived photometrical model to- 
gether with metallicities is used in chemical evolution mod- 
els to calculate the M/L-ratios of components. In the sec- 
ond part of the study (the present paper), we construct a 
consistent mass distribution model of M31. Calculating the 
rotation velocities and velocity dispersions of visible matter 
with the help of the dynamical model we found the amount 
of DM which must be added to reach an agreement with the 
observed rotation and dispersion data. 

We conclude that within uncertainties the Moore, the 
Burkert, the NFW and the Einasto distributions fit with 
observational limits at nearly all distances (Fig. 11). How- 
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ever, it is also seen that in central regions, the Moore and 
the Burkert profiles lie rather close to the upper and lower 
limits of uncertainties. The Burkert law does not fit within 
the lower limit for the mass of M 31, set by Brunthaler et al. 
(2007), cither. The NFW and the Einasto density distribu- 
tions give the best fit. 

ACDM hierarchical clustering theory is in agreement 
with dark matter density distribution in a well studied disc 
galaxy M31. Predictions of this theory also agree with ob- 
served decrease of disc sizes with rodshift (Bouwcns & Silk 
2002; Ferguson et al. 2004; Tamm & Tenjes 2006). 
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